!-------------------------------------------------------------LICENSE--------------------------------------------------------------!
!                                                                                                                                  !
!The MAP code is written in Fortran language for magnetohydrodynamics (MHD) calculation with the adaptive mesh refinement (AMR)    !
!and Message Passing Interface (MPI) parallelization.                                                                              !
!                                                                                                                                  !
!Copyright (C) 2012                                                                                                                !
!Ronglin Jiang                                                                                                                     !
!rljiang@ssc.net.cn                                                                                                                !
!585 Guoshoujing Road. Pudong, Shanghai, P.R.C. 201203                                                                             !
!                                                                                                                                  !
!This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License         !
!as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.             !
!                                                                                                                                  !
!This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of    !
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.                        !
!                                                                                                                                  !
!You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software     !
!Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.                                                   !
!                                                                                                                                  !
!-------------------------------------------------------------LICENSE--------------------------------------------------------------!

!==================================================================================================================================|
subroutine block_update(base, n_var, ro_flag, mx_flag, my_flag, mz_flag, bx_flag, by_flag, bz_flag, en_flag, po_flag)
!==================================================================================================================================|

   use parameters
   use block_define

   integer(4), intent(in) :: n_var
   logical(4), optional :: ro_flag, mx_flag, my_flag, mz_flag, bx_flag, by_flag, bz_flag, en_flag, po_flag
   type(block), pointer, intent(in) :: base

   integer(4) :: i_c, j_c, k_c, i_p, j_p, k_p, icm1, jcm1, kcm1, pos_x, pos_y, pos_z
   integer(4) :: nx_h, ny_h, nz_h, i_block, i_level, i_cpu, i_request, i_cpu_p1, i_send, i_recv
   integer(4) :: integrate_level, n_send, n_recv, n_request, nx_tmp, ny_tmp, nz_tmp, i_var
   integer(4), dimension(:), allocatable :: request
   integer(4), dimension(:, :), allocatable :: request_status
   real(8), dimension(:, :, :, :, :), allocatable :: blocks_update_var_send
   real(8), dimension(:, :, :, :, :), allocatable :: blocks_update_var_recv

   type(block), pointer :: current, parent

!----------------------------------------------------------------------------------------------------------------------------------|
!  var send flag
!----------------------------------------------------------------------------------------------------------------------------------|
   if (n_var .le. 0) then
      write(*, *)
      write(*, *) 'Warning: no variable will be updated.'
      write(*, *)
      return
   endif

!----------------------------------------------------------------------------------------------------------------------------------|
!  collect information for send and recv data
!----------------------------------------------------------------------------------------------------------------------------------|
   call mpi_allgather (n_blocks_update_send, cpu_number, mpi_integer4,                                                             &
      n_blocks_update_sendrecv, cpu_number, mpi_integer4, mpi_comm_world, error)

   n_send = sum(n_blocks_update_sendrecv(:, cpu_id_p1))
   n_recv = sum(n_blocks_update_sendrecv(cpu_id_p1, :))
   n_request = count(n_blocks_update_sendrecv(:, cpu_id_p1) .gt. 0) + count(n_blocks_update_sendrecv(cpu_id_p1, :) .gt. 0)

   allocate(request(n_request * 2), request_status(mpi_status_size, n_request * 2))
   allocate(blocks_update_var_send(n_var, nx_no_bnd, ny_no_bnd, nz_no_bnd, n_send))
   allocate(blocks_update_var_recv(n_var, nx_no_bnd, ny_no_bnd, nz_no_bnd, n_recv))

!----------------------------------------------------------------------------------------------------------------------------------|
!  collect information for sending data
!----------------------------------------------------------------------------------------------------------------------------------|
   i_send = 0
   do i_cpu = 0, cpu_number - 1
      i_cpu_p1 = i_cpu + 1
      if (i_cpu .eq. cpu_id) cycle
      n_send = n_blocks_update_sendrecv(i_cpu_p1, cpu_id_p1)
      if (n_send .eq. 0) cycle

      do i_block = 1, n_send
         i_send = i_send + 1
         i_level = 0
         current => base
         do while (blocks_update_info_send(i_level * n_dims + 1, i_block, i_cpu_p1) .ne. 0)
            pos_x = blocks_update_info_send(i_level * n_dims + 1, i_block, i_cpu_p1)
            pos_y = blocks_update_info_send(i_level * n_dims + 2, i_block, i_cpu_p1)
            pos_z = blocks_update_info_send(i_level * n_dims + 3, i_block, i_cpu_p1)
            current => current%child(pos_x, pos_y, pos_z)%pblock
            i_level = i_level + 1
            if (i_level .eq. n_levels) exit
         enddo
         do k = 1, nz_no_bnd
            k_c = k + ng
            do j = 1, ny_no_bnd
               j_c = j + ng
               do i = 1, nx_no_bnd
                  i_c = i + ng
                  i_var = 0
                  if (present (ro_flag)) then
                     i_var = i_var + 1
                     blocks_update_var_send(i_var, i, j, k, i_send) = current%ro(i_c, j_c, k_c)
                  endif
                  if (present (mx_flag)) then
                     i_var = i_var + 1
                     blocks_update_var_send(i_var, i, j, k, i_send) = current%mx(i_c, j_c, k_c)
                  endif
                  if (present (my_flag)) then
                     i_var = i_var + 1
                     blocks_update_var_send(i_var, i, j, k, i_send) = current%my(i_c, j_c, k_c)
                  endif
                  if (present (mz_flag)) then
                     i_var = i_var + 1
                     blocks_update_var_send(i_var, i, j, k, i_send) = current%mz(i_c, j_c, k_c)
                  endif
                  if (present (bx_flag)) then
                     i_var = i_var + 1
                     blocks_update_var_send(i_var, i, j, k, i_send) = current%bx(i_c, j_c, k_c)
                  endif
                  if (present (by_flag)) then
                     i_var = i_var + 1
                     blocks_update_var_send(i_var, i, j, k, i_send) = current%by(i_c, j_c, k_c)
                  endif
                  if (present (bz_flag)) then
                     i_var = i_var + 1
                     blocks_update_var_send(i_var, i, j, k, i_send) = current%bz(i_c, j_c, k_c)
                  endif
                  if (present (en_flag)) then
                     i_var = i_var + 1
                     blocks_update_var_send(i_var, i, j, k, i_send) = current%en(i_c, j_c, k_c)
                  endif
                  if (present (po_flag)) then
                     i_var = i_var + 1
                     blocks_update_var_send(i_var, i, j, k, i_send) = current%po(i_c, j_c, k_c)
                  endif
               enddo
            enddo
         enddo
      enddo
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
!  send data to the parent block which in another cpu
!----------------------------------------------------------------------------------------------------------------------------------|
   i_request = 0
   i_block = 1
   do i_cpu = 0, cpu_number - 1
      i_cpu_p1 = i_cpu + 1
      if (i_cpu .eq. cpu_id) cycle
      n_send = n_blocks_update_sendrecv(i_cpu_p1, cpu_id_p1)
      if (n_send .eq. 0) cycle
      i_request = i_request + 1
      call mpi_isend (blocks_update_info_send(:, :, i_cpu_p1), n_d_l * n_send, mpi_integer4, i_cpu, tag_int,                       &
         mpi_comm_world, request(i_request), error)
      i_request = i_request + 1
      call mpi_isend (blocks_update_var_send(:, :, :, :, i_block), n_var * nx_no_bnd * ny_no_bnd * nz_no_bnd * n_send,             &
         mpi_real8, i_cpu, tag_var, mpi_comm_world, request(i_request), error)
      i_block = i_block + n_send
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
!  receive data from the children block which in another cpu
!----------------------------------------------------------------------------------------------------------------------------------|
   i_block = 1
   do i_cpu = 0, cpu_number - 1
      i_cpu_p1 = i_cpu + 1
      if (i_cpu .eq. cpu_id) cycle
      n_recv = n_blocks_update_sendrecv(cpu_id_p1, i_cpu_p1)
      if (n_recv .eq. 0) cycle

      i_request = i_request + 1
      call mpi_irecv (blocks_update_info_recv(:, :, i_cpu_p1), n_d_l * n_recv, mpi_integer4, i_cpu, tag_int,                       &
         mpi_comm_world, request(i_request), error)
      i_request = i_request + 1
      call mpi_irecv (blocks_update_var_recv(:, :, :, :, i_block), n_var * nx_no_bnd * ny_no_bnd * nz_no_bnd * n_recv,             &
         mpi_real8, i_cpu, tag_var, mpi_comm_world, request(i_request), error)
      i_block = i_block + n_recv
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
!  waiting for mpi communication
!----------------------------------------------------------------------------------------------------------------------------------|
   call mpi_waitall (size (request), request, request_status, error)

!----------------------------------------------------------------------------------------------------------------------------------|
!  update data
!----------------------------------------------------------------------------------------------------------------------------------|
   i_recv = 0
   do i_cpu = 0, cpu_number - 1
      i_cpu_p1 = i_cpu + 1
      if (i_cpu .eq. cpu_id) cycle
      n_recv = n_blocks_update_sendrecv(cpu_id_p1, i_cpu_p1)
      if (n_recv .eq. 0) cycle

      do i_block = 1, n_recv
         i_recv = i_recv + 1
         i_level = 0
         current => base
         do while (blocks_update_info_recv(i_level * n_dims + 1, i_block, i_cpu_p1) .ne. 0)
            pos_x = blocks_update_info_recv(i_level * n_dims + 1, i_block, i_cpu_p1)
            pos_y = blocks_update_info_recv(i_level * n_dims + 2, i_block, i_cpu_p1)
            pos_z = blocks_update_info_recv(i_level * n_dims + 3, i_block, i_cpu_p1)
            current => current%child(pos_x, pos_y, pos_z)%pblock
            i_level = i_level + 1
            if (i_level .eq. n_levels) exit
         enddo

         parent => current%parent
         nx_h = nx_no_bnd / ratio
         ny_h = ny_no_bnd / ratio
         nz_h = nz_no_bnd / ratio
         nx_tmp = nx_h * (current%p_pos(1) - 1) + ng
         ny_tmp = ny_h * (current%p_pos(2) - 1) + ng
         nz_tmp = nz_h * (current%p_pos(3) - 1) + ng

         do k = 1, nz_h
            k_p = k + nz_tmp
            k_c = k * ratio
            kcm1 = (k - 1) * ratio + 1
            do j = 1, ny_h
               j_p = j + ny_tmp
               j_c = j * ratio
               jcm1 = (j - 1) * ratio + 1
               do i = 1, nx_h
                  i_p = i + nx_tmp
                  i_c = i * ratio
                  icm1 = (i - 1) * ratio + 1
                  i_var = 0
                  if (present (ro_flag)) then
                     i_var = i_var + 1
                     parent%ro(i_p, j_p, k_p) =                                                                                    &
                        sum (blocks_update_var_recv(i_var, icm1 : i_c, jcm1 : j_c, kcm1 : k_c, i_recv)) / ratio3
                  endif
                  if (present (mx_flag)) then
                     i_var = i_var + 1
                     parent%mx(i_p, j_p, k_p) =                                                                                    &
                        sum (blocks_update_var_recv(i_var, icm1 : i_c, jcm1 : j_c, kcm1 : k_c, i_recv)) / ratio3
                  endif
                  if (present (my_flag)) then
                     i_var = i_var + 1
                     parent%my(i_p, j_p, k_p) =                                                                                    &
                        sum (blocks_update_var_recv(i_var, icm1 : i_c, jcm1 : j_c, kcm1 : k_c, i_recv)) / ratio3
                  endif
                  if (present (mz_flag)) then
                     i_var = i_var + 1
                     parent%mz(i_p, j_p, k_p) =                                                                                    &
                        sum (blocks_update_var_recv(i_var, icm1 : i_c, jcm1 : j_c, kcm1 : k_c, i_recv)) / ratio3
                  endif
                  if (present (bx_flag)) then
                     i_var = i_var + 1
                     parent%bx(i_p, j_p, k_p) =                                                                                    &
                        sum (blocks_update_var_recv(i_var, icm1 : i_c, jcm1 : j_c, kcm1 : k_c, i_recv)) / ratio3
                  endif
                  if (present (by_flag)) then
                     i_var = i_var + 1
                     parent%by(i_p, j_p, k_p) =                                                                                    &
                        sum (blocks_update_var_recv(i_var, icm1 : i_c, jcm1 : j_c, kcm1 : k_c, i_recv)) / ratio3
                  endif
                  if (present (bz_flag)) then
                     i_var = i_var + 1
                     parent%bz(i_p, j_p, k_p) =                                                                                    &
                        sum (blocks_update_var_recv(i_var, icm1 : i_c, jcm1 : j_c, kcm1 : k_c, i_recv)) / ratio3
                  endif
                  if (present (en_flag)) then
                     i_var = i_var + 1
                     parent%en(i_p, j_p, k_p) =                                                                                    &
                        sum (blocks_update_var_recv(i_var, icm1 : i_c, jcm1 : j_c, kcm1 : k_c, i_recv)) / ratio3
                  endif
                  if (present (po_flag)) then
                     i_var = i_var + 1
                     parent%po(i_p, j_p, k_p) =                                                                                    &
                        sum (blocks_update_var_recv(i_var, icm1 : i_c, jcm1 : j_c, kcm1 : k_c, i_recv)) / ratio3
                  endif
               enddo
            enddo
         enddo
      enddo
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
!  update data from higher block to lower block
!----------------------------------------------------------------------------------------------------------------------------------|
   do integrate_level = max_level, 2, -1
      current => base%next
      do while (current%lv .ne. 0)
         if (current%lv .eq. integrate_level) then
            parent => current%parent
            nx_h = nx_no_bnd / ratio
            ny_h = ny_no_bnd / ratio
            nz_h = nz_no_bnd / ratio
            nx_tmp = nx_h * (current%p_pos(1) - 1) + ng
            ny_tmp = ny_h * (current%p_pos(2) - 1) + ng
            nz_tmp = nz_h * (current%p_pos(3) - 1) + ng
            do k = 1, nz_h
               k_p = k + nz_tmp
               k_c = k * ratio + ng
               kcm1 = (k - 1) * ratio + 1 + ng
               do j = 1, ny_h
                  j_p = j + ny_tmp
                  j_c = j * ratio + ng
                  jcm1 = (j - 1) * ratio + 1 + ng
                  do i = 1, nx_h
                     i_p = i + nx_tmp
                     i_c = i * ratio + ng
                     icm1 = (i - 1) * ratio + 1 + ng
                     if (present (ro_flag)) parent%ro(i_p, j_p, k_p) = sum (current%ro(icm1 : i_c, jcm1 : j_c, kcm1 : k_c)) / ratio3
                     if (present (mx_flag)) parent%mx(i_p, j_p, k_p) = sum (current%mx(icm1 : i_c, jcm1 : j_c, kcm1 : k_c)) / ratio3
                     if (present (my_flag)) parent%my(i_p, j_p, k_p) = sum (current%my(icm1 : i_c, jcm1 : j_c, kcm1 : k_c)) / ratio3
                     if (present (mz_flag)) parent%mz(i_p, j_p, k_p) = sum (current%mz(icm1 : i_c, jcm1 : j_c, kcm1 : k_c)) / ratio3
                     if (present (bx_flag)) parent%bx(i_p, j_p, k_p) = sum (current%bx(icm1 : i_c, jcm1 : j_c, kcm1 : k_c)) / ratio3
                     if (present (by_flag)) parent%by(i_p, j_p, k_p) = sum (current%by(icm1 : i_c, jcm1 : j_c, kcm1 : k_c)) / ratio3
                     if (present (bz_flag)) parent%bz(i_p, j_p, k_p) = sum (current%bz(icm1 : i_c, jcm1 : j_c, kcm1 : k_c)) / ratio3
                     if (present (en_flag)) parent%en(i_p, j_p, k_p) = sum (current%en(icm1 : i_c, jcm1 : j_c, kcm1 : k_c)) / ratio3
                     if (present (po_flag)) parent%po(i_p, j_p, k_p) = sum (current%po(icm1 : i_c, jcm1 : j_c, kcm1 : k_c)) / ratio3
                  enddo
               enddo
            enddo
         endif
         current => current%next
      enddo
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
   return
end subroutine block_update
